#include <cuda_runtime.h>
#include "freshman.h"


int recursive_reduce(int *data, int const size) {
    if (size == 1) return data[0];

    int const stride = size / 2;

    for (int i = 0; i < stride; i++) {
        data[i] += data[i + stride];
    }
    if (size % 2 != 0) {
        data[0] += data[size - 1];
    }
    return recursive_reduce(data, stride);
}

__global__ void warmup(int *g_idata, int *g_odata, unsigned int n) {
    // set tread ID
    unsigned int tid = threadIdx.x;
    // boundary check
    if (tid >= n) return;
    // convert global data pointer to the
    int *idata = g_idata + blockIdx.x * blockDim.x;
    // in-place reduction in global memory
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        if ((tid % (2 * stride)) == 0) {
            idata[tid] += idata[tid + stride];
        }
        __syncthreads();
    }
    //write result for this block to global mem
    if (tid == 0)
        g_odata[blockIdx.x] = idata[0];

}

__global__ void reduceNeighbored(int *g_idata, int *g_odata, unsigned int n) {
    // set thread id
    unsigned int tid = threadIdx.x;
    // boundary check
    if (tid >= n) return;
    // convert global data pointer to the
    int *idata = g_idata + blockIdx.x * blockDim.x;
    // in-plcae reduction in global memory
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        if ((tid % (2 * stride)) == 0) {
            idata[tid] += idata[tid + stride];
        }
        // synchronize within block
        __syncthreads();
    }
    //write result for this block to global mem
    if (tid == 0) {
        g_odata[blockIdx.x] = idata[0];
    }

}

__global__ void reduceNeighboredLess(int *g_idata, int *g_odata, unsigned int n) {
    unsigned int tid = threadIdx.x;
    unsigned idx = blockIdx.x * blockDim.x + threadIdx.x;

    int *idata = g_idata + blockIdx.x * blockDim.x;

    if (idx > n)
        return;

    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        int index = 2 * stride * tid;
        if (index < blockDim.x) {
            idata[index] += idata[index + stride];
        }
        __syncthreads();
    }
    //write result for this block to global mem
    if (tid == 0) {
        g_odata[blockIdx.x] = idata[0];
    }

}

__global__ void reduceInterLeaved(int *g_idata, int *g_odata, unsigned int n) {
    unsigned int tid = threadIdx.x;
    unsigned idx = blockIdx.x * blockDim.x + threadIdx.x;
    // convert global data pointer to the local point of this block
    int *idata = g_idata + blockIdx.x * blockDim.x;
    if (idx >= n)
        return;
    //in-place reduction in global memory
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {

        if (tid < stride) {
            idata[tid] += idata[tid + stride];
        }
        __syncthreads();
    }
    //write result for this block to global men
    if (tid == 0)
        g_odata[blockIdx.x] = idata[0];
}

int main(int argc, char **argv) {
    initDevice(0);

    bool bResult = false;

    // initialization
    int size = 1 << 24;
    printf("    with array size %d ", size);

    // execution configuration
    int blockSize = 1024;
    if (argc > 1) {
        blockSize = atoi(argv[1]);
    }
    dim3 block(blockSize, 1);
    dim3 grid((size - 1) / block.x + 1, 1);
    printf("grid %d block %d \n", grid.x, block.x);

    // allocate host memory
    size_t bytes = size * sizeof(int);
    int *idata_host = (int *) malloc(bytes);
    int *odata_host = (int *) malloc(grid.x * sizeof(int));

    // initialize the array
    initData_int(idata_host, size);

    double iStart, iElapse;
    int gpu_sum = 0;

    // device memory
    int *idata_dev = nullptr;
    int *odata_dev = nullptr;
    CHECK(cudaMalloc((void **) &idata_dev, bytes));
    CHECK(cudaMalloc((void **) &odata_dev, grid.x * sizeof(int)));

    // cpu reduction
    int cpu_sum = 0;
    iStart = cpuSecond();
//    cpu_sum = recursive_reduce(tmp, size);
    for (int i = 0; i < size; i++)
        cpu_sum += idata_host[i];

    iElapse = cpuSecond() - iStart;
    printf("cpu reduce                 elapsed %lf ms cpu_sum: %d\n", iElapse, cpu_sum);

    // kernel1: warm up
    CHECK(cudaMemcpy(idata_dev, idata_host, bytes, cudaMemcpyHostToDevice));
    CHECK(cudaDeviceSynchronize());
    // 在主机代码中调用核函数


    iStart = cpuSecond();
    warmup<<<grid, block>>>(idata_dev, odata_dev, size);
    cudaDeviceSynchronize();
    iElapse = cpuSecond() - iStart;
    cudaMemcpy(odata_host, odata_dev, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
    gpu_sum = 0;
    for (int i = 0; i < grid.x; i++)
        gpu_sum += odata_host[i];
    printf("gpu warmup                 elapsed %lf ms gpu_sum: %d<<<grid %d block %d>>>\n",
           iElapse, gpu_sum, grid.x, block.x);


    // kernel2: neighbored
    cudaMemcpy(idata_dev, idata_host, bytes, cudaMemcpyHostToDevice);
    cudaDeviceSynchronize();
    iStart = cpuSecond();
    reduceNeighbored <<<grid, block >>>(idata_dev, odata_dev, size);
    cudaDeviceSynchronize();
    iElapse = cpuSecond() - iStart;
    cudaMemcpy(odata_host, odata_dev, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
    gpu_sum = 0;
    for (int i = 0; i < grid.x; i++)
        gpu_sum += odata_host[i];
    printf("gpu reduceNeighbored       elapsed %lf ms gpu_sum: %d<<<grid %d block %d>>>\n",
           iElapse, gpu_sum, grid.x, block.x);


    // kernel3: neighbored3
    CHECK(cudaMemcpy(idata_dev, idata_host, bytes, cudaMemcpyHostToDevice));
    CHECK(cudaDeviceSynchronize());
    iStart = cpuSecond();
    reduceNeighboredLess<<<grid, block>>>(idata_dev, odata_dev, size);
    cudaDeviceSynchronize();
    iElapse = cpuSecond() - iStart;
    cudaMemcpy(odata_host, odata_dev, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
    gpu_sum = 0;
    for (int i = 0; i < grid.x; i++) {
        gpu_sum += odata_host[i];
    }
    printf("gpu reduceNeighboredLess   elapsed %lf ms gpu_sum: %d<<<grid %d block %d>>>\n",
           iElapse, gpu_sum, grid.x, block.x);


    CHECK(cudaMemcpy(idata_dev, idata_host, bytes, cudaMemcpyHostToDevice));
    CHECK(cudaDeviceSynchronize());
    iStart = cpuSecond();
    reduceInterLeaved<<<grid, block>>>(idata_dev, odata_dev, size);
    cudaDeviceSynchronize();
    iElapse = cpuSecond() - iStart;
    cudaMemcpy(odata_host, odata_dev, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
    gpu_sum = 0;
    for (int i = 0; i < grid.x; i++) {
        gpu_sum += odata_host[i];
    }
    printf("gpu reduceInterLeaved      elapsed %lf ms gpu_sum: %d<<<grid %d block %d>>>\n",
           iElapse, gpu_sum, grid.x, block.x);


    free(idata_host);
    free(odata_host);
    CHECK(cudaFree(idata_dev));
    CHECK(cudaFree(odata_dev));

    cudaDeviceReset();

    if (gpu_sum == cpu_sum) {
        printf("TEST SUCCESS!\n");
    }
    return EXIT_SUCCESS;
}